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Abstract 

We present a rigorous result on ultra-slow diffusion by solving a Fokker-Planck 
equation, which describes anomalous transport in a three dimensional (3D) 
comb. This 3D cylindrical comb consists of a cylinder of discs threaten on a 
backbone. It is shown that the ultra-slow contaminant spreading along the 
backbone is described by the mean squared displacement (MSD) of the order of 
ln(<). This phenomenon takes place only for normal two dimensional diffusion 
inside the infinite secondary branches (discs). When the secondary branches 
have finite boundaries, the ultra-slow motion is a transient process and the 
asymptotic behavior is normal diffusion. In another example, when anomalous 
diffusion takes place in the secondary branches, a destruction of ultra-slow (log¬ 
arithmic) diffusion takes place as well. As the result, one observes “enhanced” 
subdiffusion with the MSD ~ Int, where 0 < a < 1. 
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1. Introduction 


A contaminant transport in inhomogeneous media exhibits anomalous dif¬ 
fusion. This phenomenon is well established [l|, Q, Q, Q] and reviewed (see e.g., 
0, 0, 0, S, Q)- A comb model is a simple description of anomalous diffusion. 
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which however reflects many important transport properties of inhomogeneous 
media. The comb model was introduced as a toy model for understanding 
anomalous transport in low dimensional percolation clusters QQQ . It is a 
particular example of a non-Markovian phenomenon, which is also explained in 
the framework of continuous time random walks HI Q,Q- 

Anomalous diffusion on the two dimensional {2D) comb is described by the 
2D probability distribution function (pdf) P = P{x,y,t) of finding a particle 
at time t at position y along the secondary branch that crosses the backbone at 
point X. Since the transport is inhomogeneous, the diffusion on a 2D comb-like 
structure is described by the following equation [12 1 


dtP{x, y, t) = VJ{y)dlP{x, y, t) + Vyd^Pix, y, t). (1) 

Here Dx5{y) is the diffusion coefficient in the x direction, and Vy is the dif¬ 
fusion coefficient in the y direction. The ^-function in the diffusion coefficient 
in the x direction implies that diffusion occurs along the x direction at y = 0 
only. Thus, this equation describes diffusion along the backbone (at y = 0) 
where the secondary branches (fingers) play a role of traps. The comb model 
with infinite secondary branches (fingers) describes subdiffusion with the mean 
squared displacement (MSD), or the variance cr^(t), which spreads by power law 
^ [fy, [ill [13 ■ For the finite secondary branches, subdiffusion is a transient 
process until time Tq, and after a transient time scale t > Tq the transport along 
the backbone corresponds to normal diffusion with the MSD ^ t [^ . 

It is convenient to work with dimensionless variables and parameters. This 
can be obtained by the re-scaling with relevant combinations of the comb pa¬ 
rameters [Dx\ = 
coordinates are 


cm 

sec 


and [Dy] = such that the dimensionless time and 


Dlt/Dl^t, Dxx/Dy^x, Dxv/Dy ^ y/Vo, 


( 2 ) 


where D can be considered as a dimensionless diffusion coefficient for the sec¬ 
ondary branch dynamics. 

In this paper we consider anomalous ultra-slow diffusion in a three dimen- 
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Figure 1: Schematic representation of the 3D cylindrical comb. The backbone coincides with 
axis X and the discs belong to the y — z plane. 


sional cylindrical comb, which consists of a cylinder of discs threaded on the x 
axis, as it is shown in Fig.[TJ Recently, the ultra-slow phenomena were attracted 
much attention in biological search problems with long-range memories [^, 1^ . 
In the continuous time random work, ultra-slow diffusion is known as a result of 
super-heavy-tailed distributions of waiting times, see details of a discussion in 


nn 


n 


Ref. [17|, ll8| and numerical results in Ref. [19|. These heavy-tailed distributions 
result in the ratio cr^(t)/t“ —>■ 0 at t —>■ oo, which tends to zero, in contrast with 
subdiffusion, where this ratio limits to a constant value. We show that the MSD 
exhibits a logarithmic behavior in time {x^{t)) ^ ln(t), which is a result of the 
transversal branch dynamics in the 2D space. This behavior has been discussed 


in the framework of scaling arguments for the return probability 


20 


which 


are based on the fractal dimension d of the transversal branch structure and 
the spectral dimension dg, which is defined by decay of the return probab ility 
^ [§^ 21| and {x‘^{t)) ^ for ds < 2 


0 . As it is shown in Ref. [0, 


{x'^{t)) ^ Int for ds = 2 Therefore, it is instructive to present a rigorous result 
by solving the Fokker-Planck equation in the three dimensional space within 
the cylindrical comb geometry constraint, when d = 2. 
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2. Dynamics in a cylindrical comb: an infinite comb model 


We consider a 3D cylindrical comb [7|, [20|, shown in Fig.[T] in the framework 
of a standard formulation of the comb model (ED for the 3D case. Therefore, 
the random dynamics on this structure is describe by the ZD distribution func¬ 
tion P = P{x,y, z,t), where the a;-axis corresponds to the backbone, while 
the dynamics on the two dimensional secondary branches is described by the 
y and z coordinates. The diffusion equation in the dimensionless variables and 
parameters reads 


dtP = 5{:y)5{z)dlP + D{dl + dl)P. (3) 

The natural boundary conditions are taken at infinity, where the distribution 
function and its first space derivatives vanish. The initial condition is 


P{x, y, z, t = 0) = 5{x)5{y)5{z). 


( 4 ) 


2.1. Analysis in the time domain 

The formal solution of Eq. Q can be presented in a convolution form 

P{x,y,z,t)= [ G{y,z,t-t')F{x,t')dt\ (5) 

Jo 

where G{y,z,t) describes two dimensional diffusion in the secondary branches, 
while F{x, t) is a solution along the backbone. Taking into account the cylin¬ 
drical symmetry, one obtains for the y — z surface 

G(!,,,,() = G(r,f) = (6) 


where = y'^ + z'^. To define the MSD in the x direction, one needs to find a 
coarse grained distribution Pi{x,t) by integrating ([S]) over y and z taking into 
account that the element of differential area is dydz = dOrdr. From Eq. (ED, 
one obtains 


Pi{x,t) 



t') F{x,t')dt' 


( 7 ) 
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In the Laplace space, this expression establishes a relation between Pi{s) = 
and F{s) = C[F{t)]. This relation reads 

F{x, s) = sPi{x, s). (8) 


The initial condition for the coarse grained distribution is Pi{x,t = 0) = (5(x). 
Using relation (|S]), one obtains an equation for Pi{x,t). Integrating Eq. ([3]) 
over y and z, and taking into account Eqs. dil), (© and dH) , one obtains 


sA = ^]sPi(a;,s)+ (5(a;). 

Note that the Laplace transform of exists as a principal value integral 
Eourier transforming Eq. one obtains 


(9) 

22 |. 


Pi{k,s) = 


This yields the MSD in the form 


47rII 


s(47rZI + 1]) 


(x'^it)) = £ 1 - -T^Pi{k,s) 

L ak^ J A;=o 


^+200 


2ttD 


c[t- 


( 10 ) 


( 11 ) 


Taking into account that £ t ^ =t ^ and = /* + where (7 is 

an unnecessary/unimportant constant of the indefinite integration, one obtains 


{x\t)) = 


2-kD 

1 

2-kD 


ln(t) 


C 


2t: Dt 


ln(t), as t —>■ oo . 


( 12 ) 


Therefore, for the large time dynamics, ultra-slow diffusion takes place with the 
MSD being of the order of ^ ln(t). 


2.2. Consideration in the Laplace domain 

Let us consider a relation between the temporal dynamics and the dynamics 
in the Laplace space. Performing the Laplace transform of Eq. ®, one obtains 

sP = 5{y)5{z)dlP + D{dl + dl)P + Po • (13) 

Correspondingly, Eq. ([5]) reads in the Laplace domain 

P{x, y, z, s) = G{y, z, s)F{x, s) = G(r, s)F(x, s ), (14) 
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where the cylindrical symmetry is taken into account in the last term. Taking 
into account Eqs. m and (HI, one finds the solution for G{r, s) from the 
equation 

u^G” + uG' - = 0 , (15) 


where prime means a derivative over u. This is an equation for the modified 
Bessel functions Io{u) and Kq{u), where u = r^JsjD (see e.g., 23[). The solu¬ 
tion, which satisfied the boundary condition at infinity r = oo, is the modified 
Bessel function of the second kind 


G(r, s) = A • Ko(^ry/s/Dj . (16) 

It should be stressed that the Laplace inversion of Kq (^ry/s/D^ is exactly the 
solution G in Eq. 

Therefore, A = lj2TTD that satisfies the normalization condition and the initial 
condition for G{y,z,t)- Taking into account solution (fTBI) . we establish the 
relation (|S]) by integrating Eq. da over y and z. Using a property of integration 
of the modified Bessel function: 


one obtains 


uKo{au)du = \/(j? 


(17) 


P\{x,s) = F{x,s) 


dydzG{x,y,s) 


F 

D 




Fix,s) 

(18) 


which coincides exactly with the result in Eq. m, and where we also use the 
Laplace transform of Pi{x, s) = 2TrFG{r, s)rdr in Eq. ([7]). 

Now we admit an important point of the analysis: namely Eq. m cannot 
be integrated over the y and z, because G{r,s) does not exist at r = 0 (it is 
singular). Therefore, to continue the analysis, one has to return to the time 
domain and repeat the analysis for the temporal dynamics, performed in section 
\2.1\ This situation differs cardinally from the analysis for the 2D comb (HD, 
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where the finite expression for MSD can be obtained in the Fourier-Laplace 
domain. 


3. Comb dynamics with finite discs: Transition to normal diffusion 


To consider anomalous diffusion on finite combs, we consider reflected bound¬ 
ary conditions at r = R, such that drG{r = R,s) = 0, which determines the 
absence of the probability flux in the direction normal to the boundary surface. 
In this case, solution of Eq. dTSl) is found in the form of modified Bessel function 
of the first kind Io{u). Namely, it reads 


G(r, s) = Jo [(i? - r) /h , 


(19) 


which satisfies the boundary condition 
d 


dr 


{R-r)y^ 


r—R 


= h 




( 20 ) 


while for r = 0, one obtains G(0, s) = 1. This corresponds to a standard 
construction of the solution IE 


Integration of Eq. m over the y and z yields 


sPi{x, s) = dlF{x, s) -I- 1. 


( 21 ) 


Again, the relation between Pi and F can be established by integration of Io{u) 
over the y — z surface of discs, which yields 

rR r- P 


2tt 


rln 


{R — r)yjs/D dr = 2TTR^ / u/o[a(l 


u)]du. 


Jo 


Jo 


Here we used the following variable changes u = r/R and a = Ri/s/D. Then 


using another variable change w = 1 — u, one obtains two table integrals 


25| 


d 2 ra 

— / Io{w)dw — iP 

a ./n 


pi- p 

/ wIo{aw)dw = — 2 / 2 „+i(a) -/i(a) 

'o ^ ^ ^0 


Here, we are interested in the long time dynamics, when s —^ 0 and a <C 1, 
correspondingly. In this case /„(a) ~ /r(n-|- 1), and we take into account 
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the first term with n = 0 in the sum that yields Ii{a) in the squared brackets. 
Finally, one obtains for the long time asymptotic 

G(r, s)dr Ri ttR^ . (22) 

This yields the following relation 





(23) 


Substituting relation (1^51) in Eq. (EU and performing the Laplace inversion, 
one obtains the Fokker-Planck equation for normal diffusion with the diffusion 
coefficient 

dtPi = ^dlP,. (24) 

It is worth stressing that this long time diffusion takes place only for times larger 
then a transient time t > where to = R^/D. This situation is different from 


the long-time asymptotics observed in 


17|. 


This result is generic for combs with finite secondary branches (either fin¬ 
gers in the 2D comb, or discs in the 3D comb). However, the finite boundary 
conditions for the y — z discs result in the destruction of ultra-slow diffusion 
in the x direction, as well. Mathematically, this fact follows immediately from 
the Laplace image G(r,s), which depends on the boundary conditions. The 
ultra-slow motion takes place only for normal diffusion in the infinite discs. 

Note that diffusion in the side branch discs can be anomalous as well. Does 
this ultra-slow diffusion survives in this case? The answer is it does not. We 
prove this statement in the next section. 


4. Anomalous diffusion in discs 

What happens with this ultra-slow diffusion if diffusion in the transversal 
disks is anomalous and described by a memory kernel /C(t)? The comb model 

^This result can be obtained from Eq. GSli taking into account that G(r, s) 5^ 1 for the 
small argument in the limit s —>• 0. 




m now reads 


d^P = 5{y)5{z)dlP + D{dl + dl) [ IC{t - t')P{t')dt'. (25) 

Jo 


The temporal kernel IC{t) is defined in the Laplace domain through a waiting 
time pdf 4’{t) [l4- 26| 


iC = s' 0 (s)/[l - i>{s)]. 


(26) 


Repeating the analysis in the Laplace domain of Sec. 12.21 one obtains solution 
m in the form 

Gir,s)=A-KoiBr), (27) 

where B = a/ s/DlC{a) and A is a normalization constant. Therefore integration 
(fT7|) yields 

^(x,s) = ^Pi(x,s). (28) 

As already admitted above (in Sec. l 2 . 2 |) . G(r, s) is singular function at r = 0. 
Therefore, as in Eq. m, straightforward integration of Eq. (1251) over the y 
and z coordinates can be performed only in the real time domain. However, the 
function G{r = 0,t) does exists and correspondingly jC[G(r = 0,t)](s) exists as 
well, at least as a principal value integral like in Eq. Therefore, to obtain 
equation for Pi, one has to return to the time domain consideration for Pi (r, t) 
by integrating Eq. (|2^ over y and z. To be specific, let us consider subdiffusion 
in the y — z discs, described by 'ip(s) = and correspondingly with the 

memory kernel 

^(s) = 


where t is a dimensionless characteristic time scale and 0 < a < 1. In this case 
p 2 _ gOL^a Therefore, in the limit r —0 the argument Br <C 1 and Eq. 
(HZD reads for this small argument [23| 


Ko[{bs)°‘/'^] « In-- ln( 6 s), 

7 ^ 


(29) 


where 7 is an Euler constant and b = T{r'^ j. Now, we perform the Laplace 


9 








inversion at r —>■ 0 


/ ioo c\ pioQ 

KQ[{hs)^^‘^]e^^ds A5{t)\n. - Ah~^— / , 

-ioo T ^ J —ioo 


(30) 


where 4 = t/b. The last term can be presented in a form of a table integral 22| 

/ ioo 7 pioo -| 

ln(p)e^*'’dp = -— / p~^ \n{p)e^*'’dp = -. 

-200 dtb J — ioQ tb 

Finally, one obtains for t > 0 and r —>■ 0 

G{r —>• 0, t) 


Aa 1 
~ ' t 


(31) 


This result is independent of r and therefore, the limit r = 0 is correct. For this 
approximate solution, the constant A = is taken to satisfy the limit a = 1, 
which corresponds to solution (|6l) at r = 0. 

Repeating procedures of Sec. [21 namely performing first integration over 
y and z in Eq. (I25p and then the Laplace transform over time, and Fourier 
transform over x, and taking into account the result of Eq. (I28p . one obtains a 
modihcation of Eq. (fTOl) . This reads for Pi{k, s) = Cp[Pi{x, t)] 


Pi{k,s) = 


1 


(32) 


S + VaS°^k'^C[t~^) 

where Pa = aT°‘/D. Repeating the argument for the inferring Eq. ([T2l) . one 
obtains for the MSD 


{x^it)) = 


2Va P 

r(i - a) Jo {t - t'r 


OC VaP “ ln(t) . 


(33) 


Note that we omit here a term ^ jP “, which is a slower contribution to anoma¬ 


lous diffusion than the term accounted in Eq. (1551) . Here Va = 2(D„/F(2 — a) 
is a generalized transport coefficient. Subdiffusion with the transport exponent 
1 —a is dominant. Therefore, we conclude that ultra-slow diffusion ^ ln(f) takes 
place only for a = 1 that can be realized as the result of normal diffusion in the 
infinite secondary branched discs. 


5. Conclusion 

We present a rigorous result on ultra-slow diffusion by solving the Fokker- 
Planck equation in the 3D cylindrical comb geometry. It is shown that the 
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ultra-slow motion with the MSD on the x backbone is of the order of ln(t), and 
it results from normal diffusion in the secondary branched discs of the infinite 
radius. If the transport in the secondary branches is anomalously diffusive 
(subdiffusive), the anomalous transport becomes dominant in the backbone, as 
well. As the result, ultra-slow diffusion is replaced by the anomalous transport 
with the MSD ~ t^““ln(t), which is more sophisticated than usual power law 
subdiffusion, and we call it enhanced subdiffusion. This continuous transition 
from the ultra-slow motion for a = 1 to enhanced subdiffusion with 0 < a < 1 
is due to anomalous diffusion in the secondary branched dynamics, which is 
controlled by the transport exponent a. This solution can be helpful to model 
ecological processes like in generalization of an elephant random walk model 

H- 

Another important modihcation of the model is a choice of the boundary 
conditions at finite radius of the discs, which is a realistic situation. In this 
case, the physical realization of the ultra-slow transport is restricted by the 
transient time scale t < to = ttR^/D. The important and technically specific 
point of the analysis is the singularity of G(r, s) at r = 0. The integration of 
the comb equations ([3]) and (l25|l over the y and z coordinates is performed in 
the time domain, while the relation between the coarse grained pdf Pi{x, t) and 
the backbone pdf P{x, r = Off) is established in the Laplace domain. 
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